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Abstract 


In this paper, a two-phase algorithm, namely IVNS, is proposed for solving 
nonlinear optimal control problems. In each phase of the algorithm, we use a 
variable neighborhood search (VNS), which performs a uniform distribution 
in the shaking step and the successive quadratic programming, as the local 
search step. In the first phase, VNS starts with a completely random initial 
solution of control input values. To increase the accuracy of the solution 
obtained from the phase 1, some new time nodes are added and the values 
of the new control inputs are estimated by spline interpolation. Next, in 
the second phase, VNS restarts by the solution constructed by the phase 
1. The proposed algorithm is implemented on more than 20 well-known 
benchmarks and real world problems, then the results are compared with 
some recently proposed algorithms. The numerical results show that IVNS 
can find the best solution on 84% of test problems. Also, to compare the 
IVNS with a common VNS (when the number of time nodes is same in both 
phases), a computational study is done. This study shows that IVNS needs 
less computational time with respect to common VNS, when the quality of 
solutions are not different significantly. 
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1 Introduction 


Nonlinear optimal control problems (NOCP) are dynamic optimization prob- 
lems with many applications in process systems engineering, including the 
design of trajectories for the optimal operation of batch and semi-batch re- 
actors, economic systems, plasma physics, etc. [7]. 

Providing high-quality solutions with minimum computational time is the 
main issue for solving NOCPs. The numerical methods, direct [29] or indirect 
[46], usually have two main deficiencies, including low accuracy and conver- 
gence to a poor local solution. In direct methods, the quality of solutions 
depend on discretization resolution. These methods use control parametriza- 
tion to convert continuous problems to discrete problems, so they may have 
less accuracy. However, the adaptive strategies [8, 43] can overcome these 
defects, but they may be trapped by a local optimal, yet. In the indirect 
approach, the problem using Pontryagins minimum principle (PMP) is con- 
verted to two boundary value problems (TBVP) and then it can be solved by 
numerical methods such as shooting method [29]. These methods need the 
good initial guesses that lie within the domain of convergence. Therefore, 
numerical methods are not usually suitable for solving NOCPs, especially for 
large-scale and multimodal models. 

Metaheuristics as the global optimization methods can overcome these 
problems, but they usually need more computational time, though they don’t 
really need good initial guesses and deterministic rules. Several researchers 
have used metaheuristics to solve optimal control problems. For instance, 
Michalewicz et al. [34] applied floating-point Genetic algorithms (GA) to 
solve discrete time optimal control problems, Yamashita and Shima [52] used 
the classical GAs to solve the free final time optimal control problems with 
terminal constraints. Abo-Hammour et al. [1] used continuous GA for solv- 
ing NOCPs. Recently, Sun et al. [47] proposed a hybrid improved GA, for 
solving NOCPs and applied it for chemical processes. Moreover, the other 
usages of GA for optimal control problems can be found in [44, 45]. Modares 
and Naghibi-Sistani [37], proposed a hybrid algorithm by integrating an im- 
proved Particle Swarm Optimization (PSO) with a successive quadratic pro- 
gramming (SQP), for solving NOCPs. Lopez-Cruz et al. [14], applied Differ- 
ential Evolution (DE) algorithms for solving the multimodal optimal control 
problems. Recently, Ghosh et al. [22] developed an ecologically inspired op- 
timization technique, called Invasive Weed Optimization (TWO), for solving 
optimal control problems. The other well-known metaheuristic algorithms 
which are used for solving NOCPs are Genetic Programming (GP) [30], PSO 
[3, 4], Ant Colony Optimization (ACO) [48] and DE [31, 50]. 

Based on the success of the metaheuristics for solving NOCPs mentioned 
above, we propose an algorithm that use a well-known metaheuristic namely 
VNS (variable neighbourhood search) to solve NOCPs. Also, achieving a 
global optimal solution for NOCPs is another motivation for us to use a 
VNS [35]. VNS is an intelligent and metaheuristic method for solving a set 
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of combinatorial optimization and global optimization problems which uses 
neighborhood changes and uniform distributions in search procedure. Un- 
like many other metaheuristics, it is simple and requires few parameters [32]. 
Mladenovié et al. [36] proposed a general VNS for solving continuous opti- 
mization. Moreover, VNS was used for solving several optimization problem 
[25] such as mixed integer programming [26], vertex weighted k—cardinality 
tree problem [10], and scheduling problem [13]. 


In this paper, VNS uses a uniform distribution in the shaking step and 
the SQP [39], as the local search step (similar to [37]). SQP is an iterative 
algorithm for solving NLP, which uses gradient information. Furthermore, 
SQP is used for solving NOCPs alone [6, 18]. 


For performing VNS to solve an NOCP, the time interval is uniformly 
divided by using a constant number of time nodes. Next, in each of these 
time nodes, the control variable is approximated by a scalar matrix of control 
input values. Thus, an infinite dimensional NOCP is changed to a finite 
dimensional nonlinear programming (NLP). Now, we encounter two conflict 
situations: the quality of the global solution and the needed computational 
time. In other words, when the number of time nodes is increased then 
we expect the quality of the global solution to increase but we know that 
in this situation the computational time is increased dramatically. In other 
situation, we consider less number of time nodes to reduce the computational 
but we may find a poor local solution. To conquer these problems, IVNS, 
performs VNS in two phases. In the first phase of IVNS (exploration phase), 
to decrease the computational time and to find a promising solution in the 
search space, VNS uses a less number of time nodes. Next to increase the 
quality of the solution obtained from Phase 1, the number of time nodes is 
increased. Using the obtained solution in Phase 1, the values of the new 
control inputs are estimated by spline interpolation. Next, in the second 
phase of IVNS (exploitation phase), VNS uses the solution constructed by 
the above procedure, as an initial solution. A computational study in our 
numerical experiments shows that there is a significant difference between the 
computational time of IVNS and a common VNS, that uses all time nodes 
from the beginning. 


The rest of the paper is organized as follow: in Section 2, NOCPs are 
briefly introduced. In Section 3, IVNS is described. In Section 4, we provide 
more than 20 NOCPs to examine the numerical behaviour of the proposed 
algorithm. Results are compared with some numerical and metaheuristic 
methods. A computational study is carried out in Section 5 to show the 
effect of the second phase. We conclude in Section 6. 
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2 Problem formulation 


NOCPs are formulated as optimization problems by the performance index 
as the objective function and differentiate equations as constraints that called 
dynamic optimizations. There are several types of these problems e.g. track- 
ing problem, terminal control problem and time minimization problem [29]. 
We consider nonlinear bounded continuous-time control problems in which 
a vector of control functions, u, is exerted over the planning horizon [to, t;]. 
The particular problem considered is that of finding the control input vector 
u(t) € R™ that minimizes the performance index: 


min J= ¢(a(tp), ty) +f : g(x(t), u(t), t)dt (1) 
subject to: 
a(t) = f(x(t), u(t), ¢), (2) 
c(x(t), u(t), t) = 0, (3) 
d(x(t), u(t), t) < 0, (4) 
W(x(tp), ty) = 0, (5) 
x(to) = Xo, LE [to, ty]. (6) 


where z(t) € R” denotes the state vector for the system and zo € R” is 
the initial state. The functions f : R°™™ x R > R", g: R™™xR-7> 
R, c: R™™” xR = R™, d: R™™™ x R — R™, »: R” x R > R”™ and 
o~:R” x R > R are assumed to be sufficiently smooth on appropriate open 
sets. The cost function (1) must be minimized subject to dynamic (2), control 
and state equality constraints (3), control and state inequality constraints (4), 
the initial condition (6) and the final state constraints (5). 


3 Proposed algorithm 


Here, we propose IVNS for solving NOCPs. Before providing a description 
of IVNS, we introduce VNS. 


3.1 VNS algorithm 


VNS where introduced by Mladenovié and Hansen in 1997 [35] is a meta- 
heuristic algorithm which uses neighborhood changes systemically idea, both 
in the descent to local minima and in the escape from valleys which contain 
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local minima. It explores distant neighborhoods of the current incumbent 
solution, and moves from there to a new one if and only if an improvement 
is necessary. Local search method is applied repeatedly to get in the neigh- 
borhood to local optima [36]. Here, the implemented VNS in each phase has 
the following steps: 

Initialization: The time interval is divided into N;—1 subintervals using time 
nodes to,...,¢,-1 and then control input values are computed (or selected 
randomly) as control points. This can be done by the following stages: 


1. Let t, = to + kh, where h = ok =0,1,...,N;—1, be time nodes, 


where tg and ty are the initial and final times, respectively. 


2. The corresponding control input value at each time node, t,,k = 
0,...,N, —1 is an m x 1 vector, up = (ult), wes Ur having the 
following components: 

uw) =u i+ (Uright.i — Weft,i)-Ti, t= 1,2 m (7) 
a le ft,i right,i left,i}-"as gp Agneecdss 
where r; is a random number in [0,1] with uniform distribution and 
Uleft, Uright € R™ are the lower and the upper bound vectors of control 
input values, which can be given by the problem’s definition or the 
user (e.g. see the NOCPs No. 7 and 8 in Appendix, respectively). 
u = [upz]i%o' is called control input matrix. 


Evaluation: The corresponding state matrix with the control input matrix, 
u, isannx N; matrix, © = tel 5 where xp, k = 0,1,...,N¢—1, isannx1 
vector as the (k + 1)-th column of state matrix, and can approximately be 
computed by the forth Runge-Kutta method on dynamic system (2) with the 
initial condition (6). Without loss of generality, assume m = 1 (for general 
case it can be extended easily). So, the evaluation procedure is as follows: 


1 
te = Tk-1 + Elli + ale + dls + la), k=1,2,...,M:—-1 (8) 
where 
ly h 
ly =hf (xr, Ux, tk), Ig = hf(an + >> uk + 5, te) 
l h 
ly = Af (ap + 5 tk tite),  U=hf (ae +s, un + hy te) 


where uz, = u(t,) and z, = x(t,), with initial condition x(tp) = xp. To 
approximate the performance index, the composite Simpson’s method [5], is 
used. Then, the performance index in (1), J, is approximated by J as follows: 


; [341-1 ea 


I~ 9 = 6em-atw,1)+ Sf +4 S- foiti +2 y fai + f,-1) (9) 


i=l i=0 
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where f, = f(Xx, Ux, th), k = 0,1,...N,—1. If NOCP includes equality 
or inequality constraints e.g. (3) or (4), or has final state constraints, given 
by (5), then we add some penalties to the corresponding fitness value of the 
solution. Finally, we assign I(u) to u as the fitness value as follows: 


na Nie-1 ne Ny-1 


= i+ S- Mymazx{0, di(x;,u;,t;)} +30 S- Monch (@j,U;,t;) 
l=1 j=0 h=1 j=0 
Nap 
+ 5° Mapp (an,-1, tn,-1) (10) 
p=1 
where M; = [Mj1,...,Min,|’, Mz = [Moi, ..., Mon,]7? and M3 = 
[M31,... Mon, |" are big numbers, as the penalty coefficients, for cp(.,.,.), h = 


1,2,...,me, di(.,.,.), 1=1,2,...,na, and yp(.,.), p=1,2,...,ny defined in 
(3), (4) and (5), respectively. 

The fitness value in (10), can be viewed as a nonlinear objective function 
with the decision variable as u = [uo, ui,...,Un,—1]. This cost function with 
upper and lower bounds of input signals construct a finite dimensional NLP 
problem as follows: 


min J(u) = I(uo,u1,..-,UN,—-1) 
s.t 
Uleft S Uj < Uright, j=0,1,...,.M—1 (11) 


Neighborhood: VNS uses at most kmaxz neighborhoods, V,,,..., Veemaz? 1 
which r;j,2 = 1,...,kmaz is the radii of i-th neighborhood, V;, of the control 
input matrix uw. 

Shaking: In this stage, using a uniform distribution, a random direction 
matrix d € [—1,1]™** is firstly generated and then a random solution, 4, is 


selected in the k-th neighborhood, V;, by the following equation: 
j=ut+da.(r+k-—1) (12) 


where r € [0,1] is a random number, & is the index of neighborhood and a 
is the parameter of radii. 

Local search: In this stage, SQP algorithm [9, 39] is performed on the NLP 
(11), using a? = @, constructed in (12), as the initial solution when the 
maximum number of iteration is sqpmaziter. 

SQP, is an effective and iterative algorithm for the numerical solution of 
the constrained NLP problem. This technique is based on finding a solution 
to the system of nonlinear equations that arise from the first-order necessary 
conditions for an extremum of the NLP problem. Using an initial solution 
of NLP, u*, k = 0,1,..., a sequence of solutions as u*t+! = au* + d* is 
constructed, which d* is the optimal solution of the constructed quadratic 
programming (QP) that approximates NLP in the iteration k based on a*, 
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as the search direction in the line search procedure. For the NLP (11), the 
principal idea is the formulation of a QP subproblem based on a quadratic 
approximation of the Lagrangian function as L(u, A) = I(u) + ATh(u), where 
the vector A is Lagrangian multiplier and h(u) return the vector of, inequality 
constraints evaluated at u. The QP is obtained by linearizing the nonlinear 
functions as follows: 


min sat Hab) + VI(a*)Fd 
Vh(a*)?d + h(a*) <0 


Similar to [18], here a finite difference approximation is applied to compute 
the gradient of the cost function and the constraints, with the following 
components 


= 5 FSO Ajo et 13 

Ou, 5 i ; ey) 

where 6 is the double precision of machine. So, the gradient vector is 

VI = [x sides wit A Also, at each major iteration a positive definite 
= 


quasi-Newton approximation of the Hessian of the Lagrangian function, H, 
is calculated using the BFGS method [39], where \;, i = 1,...,m, is an esti- 
mated of the Lagrange multipliers. The general procedure of SQP, for NLP 
(11), is as follows: 


1. Given an initial solution u°. Let k = 0. 


2. Construct the QP subproblem (13), based on w°, using the approxima- 
tions of the gradient and the Hessian of the the Lagrangian function. 


3. Compute the new point as u*t! = a* + d*, where d* is the optimal 
solution of the current QP. 


4. Let k=k+1 and go to step 2. 


Here, in IVNS, SQP is used as the local search step, and we use the maximum 
number of iterations as the main criterion for stopping SQP. In other words, 
we terminate SQP when it converges either to local solution or the maximum 
number of SQP’s iterations is reached. 

Terminal conditions: The algorithm is terminated when the number of neigh- 
borhoods reached to kag or the difference between cost functions in two 
consecutive iterations is less than € (a given number). 

VNS algorithm is given in Algorithm 1. 
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Algorithm 1 VNS algorithm 


{Initialization} Input the number of time nodes N;, the maximum num- 
ber of iteration for SQP, sqpmaziter, a maximum number of neighborhood, 
kmax, the parameter of radii, a defined in (12), the lower and the upper 
bound vectors of control input values uje ft, Uright, an initial solution, u*, 
and «. Let k = 1. 

{Evaluation} Evaluate the fitness of the initial solution, u* and let I* = 
I(u*), where I(.) is defined in (10). 

repeat 

{Shaking} Using (12), select wu in k-th neighborhood of u*. 

{Local search} Perform SQP algorithm on the NLP (11), using u as 
the initial solution when the maximum number of iteration is sqpmaziter. 
Let u be the obtained solution, [ = I(w) and e = |I — I*|. 

if [ < I* then 

Let u* = @, I* =I and k=1. 
else 
Letk=k+1 

end if 
until k > kmaz or e <€ 
Return u* as the approximate solution, x* as the corresponding state and 
the corresponding fitness I*. 


3.2 IVNS 


We now give a new algorithm, IVNS, which is a two-phase direct metaheuris- 
tic approach. The main idea of IVNS is to find promising solution of the 
search space using the computational time as few as possible. 

IVNS has two main phases (as discussed in Section 1). In the first phase, 
we perform VNS (Algorithm 1) with a completely random initial solution 
constructed by (7). Since the main goal in this phase is to find the promising 
solution in the search space, we use a few number of time nodes. 

Next, to maintain the property of the solution given in Phase 1 and to 
increase the accurately of this solution, we add some additional time nodes. 
Thus, we increase time nodes from N;, in the Phase 1 to N;, in the Phase 
2. To use the information of the obtained solution from Phase 1 in the 
construction of the initial solution for Phase 2, we use Spline interpolation to 
estimate the values of the control inputs based on the curve obtained from 
the Phase 1. In the second phase, VNS restarts with this solution. Finally, 
IVNS is given in Algorithm 2. 


Remark 3.1. As we know, there are no general theorems on convergence of 
metaheuristics algorithm exist [28, 38]. Also, a specific theory on convergence 
of VNS does not exist, but a simple framework for global convergence of VNS 
based on attraction probabilities concept, can be found in [11]. However, we 
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Algorithm 2 IVNS 


Initialization Input ujef_ and Uright- 

{Phase 1} Perform VNS (Algorithm 1) with a random initial solution and 
using the parameters N;,, sqpmaxiter, kmax,a@ and €. (see Algorithm 1) 
{Constructing an initial solution for Phase 2} Increase time nodes 
uniformly to N;, and estimate the corresponding control input values by 
using Spline interpolation on the obtained solution from Phase 1. 

{Phase 2} Restart VNS (Algorithm 1) with the constructed initial solution 
and using N;,, s¢gpmaxiter, kmax, @ and e. (see Algorithm 1) 


mentioned that all metaheuristics are practical algorithms that are interesting 
for their numerical behaviour, [16]. 


4 Numerical experiments 


In this Section, to investigate the efficiency of IVNS, more than 20 well- 
known and real world NOCPs, as benchmark problems, are considered. These 
problems are selected with single control signal and multi control signals. 
The numerical behaviour of the algorithms can be studied from two points 
of view: the performance index and the final state constraints. Let J be 
the value of the performance index and w = [q1,..., Wel 4 defined in (5), 
and df = ||w||2 be the vector of final state constraints and the error of 7%, 
respectively. Now, the absolute errors for J and @y, are defined as follows: 


Ey=|J-—J*|, Ey =|o7 — $7 (14) 


where J* and ¢} = ||~*||2 are the best obtained solutions among the methods, 
or the exact solutions (when exist). To control the accuracy study, we now 
define a new criterion, called factor, to compare the numerical behaviour of 
the algorithms as follows: 

Ky =Ej+ Ey (15) 


where £7 and Ey, are defined in (14). Note that Ky, shows the summation 
of two important errors. Thus, based on Ky we can study the behaviour of 
algorithms on the quality and feasibility of given solutions, simultaneously. 
To solve any NOCP described in the Appendix, we must know IVNS’s 
parameters including N;:,, Ni, kmax, @, € and sqpmaziter (see Algorithm 
1), and the problem’s parameters including wie ft, Urignt and M;, i = 1, 2,3, in 
(10). To estimate the best value of these parameters, for each problem, we run 
the proposed algorithm with different values of the parameters and then select 
the best. In all NOCPs, we consider the parameters sqpmaxiter = 30, a = 
107° and kmaz = 10. The other parameters are given in the associated 
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subsection or in Table 2. Because of the stochastic nature of the proposed 
algorithm, 12 different runs were done, for each NOCP, and the best result 
are reported in Table 1. The best value of each column is highlighted in the 
bold. The reported numerical results for each algorithm included the value 
of performance index, J, the absolute error of J and E;, are defined in (14). 
The final state constraints, w = [W1,..., acl 4 the two-norm or error of the 
final state constraints, @;, the absolute error of gf and Hy, are defined in 
(14), and the factor Ky, is defined in (15). 


The algorithm was implemented in Matlab R2011la environment on a 
Notebook with Windows 7 Ultimate, CPU 2.53 GHz and 4.00 GB RAM. 
Also, to implement SQP in the proposed algorithm, as the local search, we 
used ‘fmincon’ in Matlab when the ‘Algorithm’ was set to ‘SQP’. 


In Subsection 4.1, the numerical results of IVNS are compared with exact 
solutions. Also, for comparing IVNS with metaheuristics and numerical algo- 
rithms in two Subsections 4.2 and 4.3, we consider 22 NOCPs. Their models 
are described in the Appendix, which are presented in terms of equations 
(1)-(6). The numerical results are summarized in Table 1. Details of these 
comparisons are given in the following subsection. 


4.1 Comparison with the exact solution 


Consider the nonlinear system state equations [24] 


t= 9,5 


to =u 


The cost functional to be minimized, starting from the initial states 7, (0) = 0 
and %2(0) = 1, is 


0 


The exact trajectories of the problem, from PMP, are 2j{(t) = 2 - sar 
and x3(t) = a? with the exact control signal u*(t) = teenie Also the 
exact value of the performance index is J* = 3.35. For the proposed algo- 
rithm, IVNS’s parameters are set as N;, = 15, Ni, = 21, ¢ = 10~® and the 
problem’s parameters are set as weft = —1 and Urignt = —i. The IVNS’s 


solution for the problem is J = 3.3418, thus, Ey = Ky = 0.0082. 


Figure 1 shows the graphs of the exact and the obtained trajectories, for 
x, and x2, and Figure 2 shows the graphs of the exact and the obtained 
control signals. 
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Figure 1: The exact and the obtained trajectories of (a) x; and (b) x2, for the NOCP in 
subsection 4.1 
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Figure 2: The exact and the obtained control signals for the NOCP in subsection 4.1 


4.2 Comparison with metaheuristic algorithms 


Here, six NOCPs are considered, NOCPs No. 1-6 in Appendix. The numeri- 
cal results for the first NOCP is compared with hybrid improved GA, HIGA, 
proposed in [47]. The NOCPs No. 2-4, in the Appendix are compared with 
a metaheuristic, continuous GA and CGA, proposed in [1], which gave bet- 
ter solutions than shooting method and gradient algorithm (as the indirect 
methods) [29, 12], and SUMT (as the direct methods) [18]. For NOCPs No. 
5 and 6 the results are compared with another metaheuristic, called IPSO, 
proposed in [37]. It has been shown that, for these NOCPs, IPSO was more 
accurate than some metaheuristic algorithms such as GA [42], DE [14], PSO 
[27] and some numerical methods [21, 23]. 


24 R. Ghanbari, A. Heydari and S. Nezhadhosein 


TCCR problem [47] 


The first NOCP in the Appendix is a chemical process of Temperature Con- 
trol for Consecutive Reaction, TCCR, which is an unconstrained two-state 
variable mathematical system. The objective is to obtain the optimal tem- 
perature profile that maximizes the yield of the temperature product B at 
the end of operation in a batch reactor, where the reaction A > B > C 
is occurred. The state variables, 7; and x2 are the concentration of A and 
B, respectively, and the control variable u is the temperature. The problem 
solved by HIGA [47], which was more accurate than ACO [40] and iterative 
ACO [53]. From Table 1, we can see that the numerical behaviour of IVNS 
is better than HIGA. 


VDP problem [1, 17] 


The second NOCP in the Appendix is Van Der Pol, VDP, problem which has 
two state variables and one control variable. VDP problem has a final state 
constraint, which is ~% = x1 (ty)—22(t¢)+1=0. The problem solved by CGA 
[1] and IVNS. From [1], the norm of final state constraint for the CGA equals 
ob; = 2.67 x 10~'", however, this value for IVNS equals ¢¢ = 3.04 x 107°. 
So, the factor Ky, for these methods can be seen in the sixth column of the 
Table 1. Note that the Ky, of IVNS, 3.01 x 107%, is less than CGA’s Ky, 
3.0 x 10-4. From Table 1, it is seem that IVNS can achieved more suitable 
solution than CGA. 


CRP problem [1, 29] 


The third NOCP in the Appendix is a mathematical model of Chemical Re- 
actor Problem, CRP, which has two state variables and one control variable. 
The control variable is the flow of a coolant through a coil inserted in the 
reactor that controls the first-order irreversible exothermic reaction taking 
place in the reactor. The state variables, x; and 22, are the deviations from 
the steady-state temperature and concentration, respectively. The numerical 
results of IVNS and CGA are shown in the third row of Table 1. CRP prob- 
lem has two final state constraints, w = [x1,22]". From [1], the norm of final 
state constraints for CGA, equals by = 7.57 x 10-19, when IVNS’s norm of 
final state constraints is d¢ = 2.50 x 107°. But, the corresponding Ky, of two 
methods shows that IVNS could achieve more accurate solutions than CGA. 
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FFRP problem [1, 18] 


The fourth NOCP in the Appendix is Free Floating Robot Problem, FFRP, 
which has six state variables and four control variables. It was solved by 
CGA [1]. FFRP problem has six final state constraints, wy = [7, — 4, 72, 73 — 
4,24,%5,26]’. The norm of final state constraints for IVNS is oF = 4.61 x 
10~*, however, this value, from [1], for CGA is dy = 4.65 x 107%. From Table 
1, we can see the numerical behaviour of IVNS is better than CGA, also it 
is clear that the obtained values of J, Ez, of, Hy and Ky from IVNS are 
better than CGA. 


CSTCR problem [37] 


The fifth NOCP in the Appendix is a model of a nonlinear Continuous Stirred- 
tank Chemical Reactor, CSTCR. It has two state variables 71(t) and x2(t), as 
the deviation from the steady-state temperature and concentration, and one 
control variable u(t), which represents the effect of the flow rate of cooling 
fluid on chemical reactor. The objective is to maintain the temperature and 
concentration close to steady-state values without expending large amount 
of control effort. Also, this is a benchmark problem in the handbook of 
test problems in local and global optimization [20], which is a multimodal 
optimal control problem [2]. It involves two different local minima. The 
values of the performance indices, for these solutions, equal 0.244 and 0.133. 
The numerical results of IVNS, with the parameters in Table 2, are compared 
with IPSO [37], and numerical methods in [2, 14]. From the results of the 
fifth row of Table 1, we can see that IVNS is the best. 


MSNIC problem [37] 


In the sixth NOCP in the Appendix, a Mathematical System with Nonlinear 
Inequality Constraint, MSNIC, is considered. It includes an inequality con- 
straint, d(x,t) = r2(t) + 0.5 — 8(t — 0.5)? < 0. From the sixth row of Table 
1, we can see that the obtained value of the performance index, for IVNS 
is J* = 0.1720, which is better than IPSO’s, 0.1727, and other numerical 
methods given in [23, 33]. 


4.3 Comparison with numerical algorithms 


In this subsection, for NOCPs no. 7-22, the results of IVNS are compared 
with some numerical methods such that SQP [18], SUMT [18], Bézier [21], 
HPM [15], DTM [41] and ADM [19]. Usually, for these methods the final 


26 R. Ghanbari, A. Heydari and S. Nezhadhosein 


state constraints are not reported. But these values are reported for IVNS 
in Table 1. 


Comparison with Bézier [21] 


The NOCP No. 7, in the Appendix, has exact solution, i.e. the exact value of 
performance index equals J* = —5.5285 [49]. This problem has an inequality 
constraint as d(x,t) = —6 — 21(t) < 0. It has been solved by a numerical 
method, proposed in [21], called Bézier, and the proposed algorithm, IVNS, 
with the parameters in Table 2. From seventh row of Table 1, the obtained 
value of the performance index from IVNS is better and more accurate than 
Bézier method. 


Comparison with HPM [15], DTM [41] and ADM [19] 


In this subsection, the results of [VNS with the parameters given in Table 2, 
are compared with HPM [15], DTM [41] and ADM [19]. For NOCP No. 8 in 
the Appendix, which is a constraint nonlinear model, the numerical results 
are compared with HPM. This NOCP has a final state constraint as 


wy=x-0.5=0. 


From [15], the norm of final state constraint for HPM is dy = 4.2 x 10~°, 
however, this value for IVNS equals ¢% = 6.83 x 107'*. From Table 1, it is 
clear that the obtained values of the performance index, the norm of final 
state constraint and Ky from IVNS are better than HPM’s. 

The problem No. 9 in the Appendix is a linear quadratic optimal control 
which has been solved by two numerical methods, DTM [41] and ADM [19]. 
Using the approximate values of k(t), which is used to achieve the optimal 
control signal by linear feedback control as u(t) = —k(t)a(t), the performance 
index could be calculated. The exact solution, from PMP, equals J* = 0.1929. 
From Table 1, the values of Ey and Ky, for IVNS, with the same number of 
points, N,, = 15, equals 0.0052, which is less than DTM and ADM methods, 
(0.0087). 


Comparison with SQP and SUMT 


For NOCPs No. 10-22 in the Appendix, the numerical results of IVNS (the 
parameters are given in Table 2) are compared with SQP and SUMT meth- 
ods. All these problems are described in [18]. For SQP and SUMT, the status 
of the final state constraints were not reported, so, we replaced the values of 
gy instead of Ey, in Table 1. Also, in computation of the factor, Ky, the 
values of Fy, for SQP and SUMT methods are considered to be zero. The 
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results (given in Table 1) show that IVNS could find more accurate results 
for performance index J, and the factor Ky, perspective. 


Table 1: The best of numerical results for 12 different runs of NOCPs described in 
Appendix 


Problem Algorithm 7 By Ey Ky 
TCOCR HIGA[47] 0.61046 2.0 x 1075 = 2.0 x 1075 

IVNS 0.61048 0 — 0 
VDP CGA [1] 1.7404 3.0 x 10-7 i) 3.0 x 10-7 

IVNS 1.7401 () 3.01 x 1079 3.01 x 1079 
GRP CGA[I] 0.0163 4.0 x 10-4 0 4.0 x 10-7 

IVNS 0.0159 te) 2.42 x 1078 2.42 x 1078 
FFRP CGATI] 83.63 1772 0.0042 17.7242 

IVNS 65.91 0 0 0 
CSTCR IPSO [37 0.1354 0.0024 = 0.0024 

[2] J € (0.135, 0.245] 0.0020 = 0.0020 

[14] J € [0.1358,0.1449] 0.0028 = 0.0028 

IVNS 0.1328 2.0 x 1074 = 2.0 x 1074 
MSNIC TPSO [37 0.1727 0.0007 = 0.0007 

[23] 0.1816 0.0096 — 0.0096 

[33] 0.1769 0.0049 = 0.0049 

IVNS 0.1720 0 = oO 
NOCP no. 7 Bézier [21] =5. 3898 0.1387 = 0.1387 

IVNS —5.5082 0.0203 = 0.0203 
NOCP no. 8 HPM [15 0.2353 0.0338 4.20 x 10-6 0.0338 

IVNS 0.2015 0 0 0 
NOCP no. 9 DTM [41 0.2016 0.0087 = 0.0087 

ADM [19 0.2016 0.0087 — 0.0087 

IVNS 0.1877 0.0052 = 0.0052 
NOCP no. 108 SUMT [18 5.15 x 10-6 5.14 x 10-6 = 5.14 x 10-6 

SQP [18 6.57 x 10~8 6.56 x 1078 6.56 x 1078 

IVNS 6.57 x 10711 i) = te) 
NOCP no. 11° SUMT [i8 1.7980 0.0791 = 0.0791 

SQP [18 1.7950 0.0761 = 0.0761 

IVNS 1.7189 0 = 0 
NOCP no. 12° SUMT [i8 0.1703 0.0223 = 0.0223 

SQP [18 0.2163 0.0683 — 0.0683 

IVNS 0.1480 0 — 0 
NOCP no. 13° SUMT [i8 3.2500 0.3507 NR® 0.3507 

SQP [18 3.2500 0.3507 NR 0.3507 

IVNS 2.8993 te) 7.49 x 10710 7.49 x 10710 
NOCP no. 14° SUMT [i8 —0.2490 0.001 NR 0.001 

SQP [18 —0.2490 0.001 NR 0.001 

IVNS —0.2500 ts) 2.6 x 10710 2.6 x 10719 
NOCGP no. 15° SUMT [18 0.0167 6.0 x 10-4 NR 6.0 x 10-4 

SQP [18 0.0168 7.0 x 1074 NR 7.0 x 10-4 

IVNS 0.0161 fe) 3.42 x 1079 3.42 x 1079 
NOCP no. 16 SUMT [18 3.7700 0.4648 NR 0.4648 

SQP [18 3.7220 0.4168 NR 0.4168 

IVNS 3.3052 te) 3.35 x 1078 3.35 x 1078 
NOCP no. 172 SUMT [18 9.29 x 10-4 3.0 x 10-8 NR 3.0 x 10-5 

SQP [18 1.01 x 1073 8.4 x 1075 NR 8.4 x 1075 

IVNS 9.26 x 1074 ts) 6.66 x 10719 6.66 x 10719 
NOCP no. 18° SUMT [18 2.2080 0.2079 NR 0.2079 

SQP [18 2.2120 0.2119 NR 0.2119 

IVNS 2.0001 te) sal x io" 5.01 x 10711 
NOCP no. 19° SUMT [18 —=8.8690 0.0002 NR 0.0002 

SQP [18 —8.8690 0.0002 NR 0.0002 

IVNS —8.8692 te) 6.89 x 1079 6.89 x 1079 
NOCP no. 20° SUMT [18 0.0368 0.0042 = 0.0042 

SQP [18 0.0368 0.0042 = 0.0042 

IVNS 0.0326 0 — oO 
NOGP no. 21° SUMT [18 76.83 12.11 NR 12.11 

SQP [18 77.52 12.80 NR 12.80 

IVNS 64.72 () 1.46 x 1074 1.46 x 1074 
NOCP no. 22° SUMT [18 0.3428 0.0670 NR 0.0670 

SQP [18 0.3439 0.0681 NR 0.0681 

IVNS 0.2758 te) 0.0021 0.0021 


™Not Reported. 
> We here consider, Ey, = ¢f for IVNS, and for SQP and SUMT methods, Ey =0 
(since the values were not reported, we consider the best possible situation for SQP and SUMT). 


Table 1 shows that IVNS was 100 percent successful in point of view the 
performance index, numerically. The associated values of EF; for IVNS are 
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zero for all test problems. It shows that IVNS provides robust results with 
respect to the other methods. 

To have a more careful comparison, we computed the Gap between the 
performance index’s value of the algorithms and the best obtained perfor- 
mance index’s value. In other words, let J be the obtained value of the 
performance index of an algorithm. Now, similar to [51], we define the Gap 
as follows: ys 

Gap(J) =| 
From Table 1, the mean values of Gap for IVNS, SQP and SUMT, on NOCPs 
No. 10-22, are 0, 7.69e + 3 and 6.02e + 3, respectively. Thus it is obvious 
that, IVNS gave more better solution in comparison with SQP and SUMT. 
We believe that this is due to the fact that IVNS tries to find the global 
solution but SQP and SUMT didn’t escape from a local minimum. 

To compare with the CGA (as a global search algorithm), from Table 
1, we see that the mean values of the Gap for CGA is 0.0981. Thus, we 
can see IVNS is 100 percent better than CGA from Gap perspective. This 
result shows that IVNS’s estimations of global minimal is better than CGA’s 
estimation. Therefore, based on these numerical study, we can conclude that 
IVNS outperforms than CGA. 

The mean values of violation of the norm of the final state constraints, of, 
for IVNS is 1.16 x 107-4. Therefore, it is evident that IVNS is more robust. 
Also, the mean value of ¢f for IVNS and CGA are 1.53x 1074 and 1.55 107%, 
respectively, on NOCPs no. 2-4. Thus, we can say that the feasibility of the 
solutions given by IVNS and CGA are competitive. Therefore, it is seen that 
IVNS could provide very suitable solutions with respect to the optimality 
and feasibility criteria. Also, the mean of the factor, Ky, for IVNS equals 
1.28 x 10-3. For NOCPs No. 10-22 the mean of factor for IVNS, SQP 
and SUMT equals 1.76 x 10-4, 1.0768 and 1.0272, respectively. Therefore, 
we can say that IVNS outperform well-known numerical methods. Since, 
the computational times of the most algorithms were not reported thus we 
didn’t give the computational times of IVNS in Table 1. But, the details of 
the computational time of IVNS is given in Table 3 that will be discussed in 
Section 5. 


| (16) 


5 Comparison with a common VNS 


The main idea for proposing a two-phase algorithm is to decrease the required 
computational time in solving NOCPs. So, we focus on investigating of the 
influence of the second phase in IVNS. To compare the IVNS with a common 
VNS, the number of time nodes are selected same in both phases. In common 
VNS, only the first phase of IVNS, which the number of time node equal 
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Table 2: The parameters of IVNS for NOCPs described in the Appendix 


Problem Ule ft Uright Nz, Ni, E M; 
TCCR 298 398 411 15 10° = 
VDP -0.5 2 31 151 10-6 7 

CRP 5 2 of. bk 0 (1, 1)” 
FFRP 215 10 SL. 6 105? [es TOG 4 
CSTCR 0 5 31 5 107% = 
MSNIC  -20 20 a br. “io i 

no. 7 -2 2 21 131 10-7 1 

no. 8 -2 2 31 91 107° 1 

no. 9 me 3 11 15 10-8 _ 

no. 10 -3 3 21 51 10-6 = 

no. 11 -2 2 31 91 10° 1 

no. 12 -20 20 315i te- 1 

no. 13 -4 3 31 75 10-6 100, 100)” 
no. 14 il 1 St. 7 “10-8 1000 

no. 15 =) 2 7 6410S 10-8 100, 100] 
no. 16 —7 1 a BL 10-* 100, 100)” 
no. 17 -1 1 fi 6 -35~=C 10-8 [10, 10]” 
no. 18 -5 5 31 151 10-6 (10, 10]7 
no. 19 -30 30 31° AW 10-* 100, 100)" 
no. 20 -1 1 31 171 107° — 

no. 21 215 10 at TO 0 cn, TO 
no. 22 -15 10 a. Ot 0" (10... 1, 


N:,, is applied. For these methods, 35 different runs, for each NOCP in 
the Appendix, were made with the same parameters. The influence of these 
methods investigated for these NOCPs on the dependent outputs consist of 
performance index, J, the factor, df and required computational time, Time. 
The results are given in Table 3. 

From Table 3, we observe that the two-phase method has no significant 
effect on J, of. But the two-phase method, IVNS, needs less computational 
time than the common VNS, significantly (except NOCP No. 16). Therefore, 
based on this computational study, we can conclude that the usage of two- 
phase VNS can decrease the computational time, significantly, without loss 
of quality of solution. 


6 Conclusion 


In this paper, a two-phase algorithm, namely IVNS, was proposed for solving 
NOCPs. In each phase of the algorithm, we used a VNS, which performed 
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Table 3: The best numerical results for NOCPs in Appendix, using IVNS and common 
VNS 


IVNS VNS 
Problem J of Time J of Time 
TCCR 0.6105 — 4.1496 0.6107 = 4.2482 
VDP 1.7401 3.04x 10-9 375.69 1.7513 1.42x10-® 413.28 
CRP 0.0159 2.50x 10-8 = 78.09 0.0164 3.12x 10-9 = 112.05 
FFRP 65.91 4.61x 10-4 264.62 50.31 8.17 x 1073 —- 285.13 
CSTCR 0.1328 — 48.82 0.1116 — 52.83 
MSNIC 0.1720 — 10.49 0.1725 — 29.82 
no. 7 —5.5082 — 42.27 —5.5012 — 65.81 
no. 8 0.2015 6.83 x 1071! 11.18 0.2012 4.21 x 10719 12.24 
no. 9 0.1877 — 3.1278 0.1899 — 5.6636 
no. 10 6.57 x 10-14 — 3.7440 2.88 x 10-14 — 3.9624 
no. 11 1.7189 = 119.94 1.7152 — 139.55 
no. 12 0.1480 — 41.38 0.1486 — 54.35 
no. 13 2.8993 7.49 x 10-19 39.04 2.8935 3.41x 10-9 38.36 
no. 14 —0.2500 2.60 x 10-19 52.61 —0.2498 1.52x 10-8 = 93.10 
no. 15 0.0161 3.42x 10-9 124.02 0.0162 4.03 x 107!9 154.65 
no. 16 3.3052 3.35 x 10-8 — 137.85 3.3051 1.03 x 10-19 111.07 
no. 17 9.26x10-4 6.66 10-19 144.16 9.81x10-4 8.35x10-8 178.23 
no. 18 2.0001 5.01 x 10-1! 35.10 2.0001 2.13 x 10-12 120.07 
no. 19 —8.8692 6.89 x 10-9 114.30 —8.8692 7.13x 10-9 129.02 
no. 20 0.0326 — 42.69 0.0326 — 64.23 
no. 21 64.72 1.46x 10-4 =: 145.68 56.54 4.74x 107% = 148.01 
no. 22 0.2758 0.0021 135.25 0.2765 0.0038 217.65 


a uniform distribution in the shaking step and the SQP, as the local search 
step. In the first phase, VNS started with a completely random initial solu- 
tion of control input values. To increase the accuracy of the solution obtained 
from Phase 1, the some new time nodes were added and the values of the new 
control inputs were estimated by Spline interpolation. Next, in the second 
phase, VNS restarted by the solution constructed by Phase 1. Finally, we im- 
plemented the proposed algorithm on more than 20 well-known benchmarks 
and real world problems, then the results were compared with some recently 
proposed algorithms. The numerical results showed that IVNS could found 
mostly better solution than other proposed algorithms. Also, to compare of 
IVNS with a common VNS a computational study was done that showed 
that IVNS needed less computational time with respect to a common VNS. 
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pendix 


following NOCPs are described using eqns (1)-(6). 


[47, 53, 40] (TCCR) ¢ = 22, to =0, ty = 1, f = [-4000exp(—2500/u) 
x2, 4000exp(—2500/u)x? — 620000exp(—5000/u)ae]7, d = [298 —u, u— 
398)", xo = (1, 0]7. 


" (1, 17] (VDP) g= 5 (a7 + 25 + u?), to = 0, ty = 5, f — [t2, —22 + (1 = 
x7)ax9 =F ul”, xo = [1,0]7, p =2£,—-%2+ 1. 


. [1, 29] (CRP) g = $(a7+23+0.1u?), to = 0,t¢ = 0.78, f = [v1 -2(a1 + 
0.25) + (xq + 0.5)erp(2521/(x1 + 2)) — (a1 + 0.25)u, 0.5 — rg — (a2 + 
0.5)exp(2521 /(x1 + 2))|", zo = [0.05, 0)", » = [21, xe]”. 


. [l, 18] (FFRP) g = $(u? + u3 + v3 + u}),to = Oty = 5,f = 
[v2, ((ui + ug) cosas — (ug + u4) sinzs5)/M, x4, ((ui + ug) sins + (ug + 
ua) cos x5) /M, Tv, (D(ur +u3)—Le(us +ua))/I]", 20 = [0, 0, 0, 0, 0, we 
wb = [a1 — 4, 22,23 — 4,24,25, 26), M = 10,D =5,1 = 12, Lh, =5. 


9 


. [37] (CSTCR) g = 2? +.23+0.1u?, to = 0,t7 = 0.78, f = [-(2+u)(214 
0.25) + (#2 +0.5)eap(25a1/(a1+2)),0.5—a2 — (42g +0.5)exp(2521/(a1+ 
2))|7, 29 = (0.09, 0.09)". 


. [37] (MSNIC) @ = 23,to = 0,tp 1,f [v2,-%2 + u,v? + 22 + 
0.005u2]7, d = [-(20—u)(20-+4), 22 +0.5—8(t—0.5)2]?, v9 = [0,-1, 0]. 


. [21] g = 2x1, to = 0,ty = 3, f = [x2,u]?,d = [-(2— u)(2 + u), -6 
£1)" , 20 = (2, 0]”. 


36 


15. 


16. 


I. 


18. 


19. 


20. 


21. 


22. 


R. Ghanbari, A. Heydari and S. Nezhadhosein 


15] g=u?,to =0,ty =1,f = 52? sing +u, x =0,p=2-0.5. 


18] g cos“ u, tp =0,ts =7, f=sin >, ro = 
18] g = 5(@7 + 23 + u?), to = 0,ty = 5, f = [e2,-a1 + (1 — aj)aa + 
ul? ,d = —(x2 + 0.25), xo = [1, 0]7. 


18] g = a7 + 23 + 0.005u?, to = 0,t7 = 1,f = [x2,-2o + ul’, d = 
—(20 — u)(20 + u), 0.5 + 22 — 8(t — 0.5)?]7, 29 = [0, —1]*. 


18] g = fu, to = 0, te = 2, f = [20, ul”, eo = [1,1]? , b = [21, x0]”. 
18] g = ae o = 0,t7 = 1,f = (22, u]?,d = —(1—u)(1+ 4), 29 = 
0,0]7,% = 


18] g= er = 0,t¢ = 0.78, f = [-2(@1 +0.25) + (a2 
0.5)eap(25a1 /(a1+2))— (a1 +0.25)u, 0.5—22— (a2 +0.5)exp(25a1 /(a1 + 
2)) Tae = (0.05, 0}”, a = [ei wal. 


[18] g = Zu?, to = 0, ty = 10, f = [cosu — x2, sinu]’,d = —(m — u)(m + 
u), £9 = [3.66, —1.86]”, = (21, xe]”. 


[18] g = (a7 + 23),to = 0,t¢ = 0.78, f = [—2(a1 + 0.25) + (x9 
Boe xp(2 are Get 2) Osher, 05 — ayn US ep nea aa 
2))]?,d=—-(1—u)(1 +4), zo = [0.05, 0]7, b = [21, z9]". 


18 co) _ £3, to 0, ty Lf [r2,u alee d= ty — 1.9, Zo = 

0,0, 0]? = [21,22 + 17. 

18] $ = —a3,to = 0,ty = 5, f = [o2,-2+ *,—0.0lu]?,d = —(30 - 

u)(30 + u), 29 = [10, —2, 10]7, b = [x1, 29]". 
1 
gu 


18] d= (ws — 1)’ + 03 + 5 g= 
sin u, sin ul’, 29 = [0, 0, 0)”. 


ww, io = 0,t¢—5, f = [ee coat, eg 


18] g = (uj +u3 +.u3 +. uJ), to = 0, ty = 5, f = [x2, (ur +us) cos a5 — 
(ug+u4) sin as) /M, va, ((uyt+ug) sin v5-+(u2+u4) cos v5) /M, x6, (D(uy+ 
us) _ Le (ua a ua))/I]", xo — [0, 0, 0,0, 0, 0]7, w = [ary _ 4,22, %3 — 


4,24, 5 — =,26]7, M =10,D =5,I1 =12, L. =5. 


[18] g = 4.5(a3+22)+0.5(uf+u3), to = 0, te = 1, f = (9x4, 9x5, 926, 9(u1 
+17.2523), 9ue, ~9(tin 27.075 6cr9-4 2arsarg) sto] ,ro = [0, 22, 0, 0, —1, 0]? 
Ww _ [v1 _ 10, x2 = 14, 3,04 — 2.5, U5, ae)". 
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